#####################################################################
## Replication file for                                           
## When Do Citizens Consider Political Parties Legitimate?  
##                                   
## Author: Ann-Kristin Kölln (ann-kristin.kolln@gu.se)                                                     
## 26th May 2023                                                 
##                
#####################################################################


# R version 4.3.0 (2023-04-21) -- "Already Tomorrow"; Platform: x86_64-apple-darwin20 (64-bit)



## load data European-wide survey Party Legitimacy
data_long2<-read.csv("PartyLeg_crosssec.csv", header=TRUE)


data_long2$Edu_o<- as.factor(data_long2$Edu_o)
data_long2 $Edu_o <- factor(data_long2 $Edu_o, levels = c("low", "medium","high"))

data_long2$Extremity_SelfLRN <- as.factor(data_long2$Extremity_SelfLRN)
data_long2 $Extremity_SelfLRN <- factor(data_long2 $Extremity_SelfLRN, levels = c("centre", "moderate left/right","left/right wing"))

data_long2$gov_exp <- as.factor(data_long2$gov_exp)
data_long2$newer_party <- as.factor(data_long2$newer_party)


## install and load required packages
install.packages(c("radiant.data", "plyr","ggplot2", "Rmisc", "tidyr", "ggpubr", "lme4", "lmerTest", "sjPlot", "effects", "psych", "ltm", "lavaan"))

library(radiant.data)
library(plyr)
library(ggplot2)
library(Rmisc)
library(tidyr)
library(ggpubr)
library(lme4)
library(lmerTest)
library(sjPlot)
library(effects)
library(psych)
library(ltm)
library(lavaan)


### Descriptives ###

### weighted means and sd of index and of individual items
weighted.mean(data_long2$Leg_01, data_long2$wFactor, na.rm=TRUE)
weighted.sd(data_long2$Leg_01, data_long2$wFactor, na.rm=TRUE)

weighted.mean(data_long2$Leg_SameRights, data_long2$wFactor, na.rm=TRUE)
weighted.sd(data_long2$Leg_SameRights, data_long2$wFactor, na.rm=TRUE)

weighted.mean(data_long2$Leg_DiscussProposals, data_long2$wFactor, na.rm=TRUE)
weighted.sd(data_long2$Leg_DiscussProposals, data_long2$wFactor, na.rm=TRUE)

weighted.mean(data_long2$Leg_AcceptDecisions, data_long2$wFactor, na.rm=TRUE)
weighted.sd(data_long2$Leg_AcceptDecisions, data_long2$wFactor, na.rm=TRUE)             
               


## producing Figure 1 by summarizing data and plotting country-specific plots that are compiled

# setting theme
theme_set(theme_sjplot())


leg_parties <- summarySE(data_long2, measurevar="Leg_01", groupvars=c("PartyName", "CountryName"), na.rm=TRUE)
leg_parties2<-na.omit(leg_parties)

leg_partiesDE<-leg_parties2[leg_parties2$CountryName == "Germany",]
leg_partiesDK<-leg_parties2[leg_parties2$CountryName == "Denmark",]
leg_partiesES<-leg_parties2[leg_parties2$CountryName == "Spain",]
leg_partiesFR<-leg_parties2[leg_parties2$CountryName == "France",]
leg_partiesNL<-leg_parties2[leg_parties2$CountryName == "Netherlands",]
leg_partiesPT<-leg_parties2[leg_parties2$CountryName == "Portugal",]
leg_partiesUK<-leg_parties2[leg_parties2$CountryName == "UK",]


p_leg_DE<- ggplot(leg_partiesDE, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,0.5,0,1), "cm"))

p_leg_DK<- ggplot(leg_partiesDK, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,0.5,0,1), "cm"))
	
p_leg_ES<- ggplot(leg_partiesES, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,1,0,1), "cm"))

p_leg_FR<- ggplot(leg_partiesFR, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,1,0,1), "cm"))

p_leg_NL<- ggplot(leg_partiesNL, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,1,0,1), "cm"))
	
p_leg_PT<- ggplot(leg_partiesPT, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,0.5,0,1), "cm"))

p_leg_UK<- ggplot(leg_partiesUK, aes(x=PartyName, y= Leg_01)) + 
    geom_errorbar(aes(ymin=Leg_01-sd, ymax= Leg_01+sd), width=.5) +
    geom_point(size=2)+
    xlab("")+
    ylab("Legitimacy Index")+ 
    ylim(0.2, 1.1) +
	theme_bw() + theme(plot.margin = unit (c(1,1,0,1), "cm"))
	
	

# all in one plot to produce Figure 1
jpeg("Figure1.jpeg", width= 10, height= 10, units='in', res=300)
ggarrange(p_leg_DK, p_leg_FR, p_leg_DE, p_leg_NL, p_leg_PT, p_leg_ES, p_leg_UK, labels = c("Denmark","France","Germany","Netherlands","Portugal","Spain","UK"), ncol = 2, nrow = 4)
dev.off()


# democratic behavior: weighted mean and sd
weighted.mean(data_long2$DemBeh_01, data_long2$wFactor, na.rm=TRUE)
weighted.sd(data_long2$DemBeh_01, data_long2$wFactor, na.rm=TRUE)



### Hypothesis tests ###

# evaluating the hierarchical model
Model0<-lmer(Leg_01~ 1 + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model0a<-lmer(Leg_01~ 1 + (1|Party), data=data_long2, weights = wFactor)

anova(Model0, Model0a)

tab_model(Model0, Model0a)

## input for Table 1

# H1: gov -> Leg
Model2<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)

# H2: inst -> Leg
Model3<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + newer_party + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)

# H3: extr -> Leg
Model4<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_LR  +  (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)


# H4: DemBeh -> Leg
Model5<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + DemBeh_01 + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)

# Full model
Model6<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + newer_party + IdeolModeration_LR + DemBeh_01 + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)

## Producing Table 1
tab_model (Model2, Model3, Model4, Model5, Model6, show.ci = FALSE, show.se = TRUE, emph.p =FALSE, digits.p = 2, file ="Table1.doc")


### postestimation analysis

allEffects(Model4)
allEffects(Model5, xlevels = list(DemBeh_01= seq(0.24,0.82, 0.02)))




#### Supplementary Information ####

## Part B ##


### evaluate indices: legitimacy and democratic behaviour

## Legitimacy

# correlations between legitimacy variables
Leg<-data_long2[c("Leg_SameRights", "Leg_DiscussProposals", "Leg_AcceptDecisions")]

cor(Leg, use = "pairwise.complete.obs")

# calculating Cronbach's alpha 
cronbach.alpha (Leg, na.rm=TRUE) 

# PCA
Leg_NA<-na.omit(Leg)

Leg_NA.pca<-prcomp(Leg_NA, center = TRUE, scale =TRUE) 
summary(Leg_NA.pca)
Leg_NA.pca$rotation 


## Democratic behaviour

# correlations between variables
DemBeh<-data_long2[c("ComplyLaws", "EndorseDemo")]
cor(DemBeh, use = "pairwise.complete.obs")

## calculating Cronbach's alpha 
cronbach.alpha (DemBeh, na.rm=TRUE)

## PCA
DemBeh_NA<-na.omit(DemBeh)
DemBeh_NA.pca<-prcomp(DemBeh_NA, center = TRUE, scale =TRUE) 
summary(DemBeh_NA.pca)
DemBeh_NA.pca$rotation 






### validate legitimacy index against other constructs

##correlations with sympaty and vote probabilities
cor.test(data_long2$Sympathy, data_long2$Leg_01)
cor.test(data_long2$VoteProbability, data_long2$Leg_01)


# test differences in means with a binary version of legitimacy
data_long2$Leg_bin[data_long2$Leg_01 < 0.51]<-0
data_long2$Leg_bin[data_long2$Leg_01 > 0.50]<-1

ddply(data_long2, ~ Leg_bin, summarise, mean = mean(VoteProbability, na.rm=TRUE), sd = sd(VoteProbability, na.rm = TRUE), N = length(VoteProbability))
fit_vote<-aov(VoteProbability ~ Leg_bin, data=data_long2)
summary(fit_vote)

table(data_long2$Leg_bin, data_long2$Sympathy)
ddply(data_long2, ~ Leg_bin, summarise, mean = mean(Sympathy, na.rm=TRUE), sd = sd(Sympathy, na.rm = TRUE), N = length(Sympathy))
fit_symp<-aov(Sympathy ~ Leg_bin, data=data_long2)
summary(fit_symp)





### Confirmatory Factor analyses

## Sympathy
leg_1<- 'leg=~ Leg_SameRights + Leg_DiscussProposals + Leg_AcceptDecisions
			eval=~Sympathy'
			
fit_leg_1<- cfa(leg_1, data=data_long2)

# summary
summary(fit_leg_1, fit.measures= TRUE, standardized=TRUE) # really good values



## Vote probabilities
leg_2<- 'leg=~ Leg_SameRights + Leg_DiscussProposals + Leg_AcceptDecisions
			eval=~ VoteProbability'
			
fit_leg_2<- cfa(leg_2, data=data_long2)

# summary
summary(fit_leg_2, fit.measures= TRUE, standardized=TRUE) # really good values again



## Table B1: producing input
TableB1 <- ddply(data_long2, c("PartyName"), summarise,
               N    = sum(!is.na(Leg_01)),
               mean = weighted.mean(Leg_01, w= wFactor, na.rm=TRUE),
               sd   = sd(Leg_01, na.rm=TRUE),
               se   = sd / sqrt(N))
               
write.csv(TableB1, file = "TableB1.csv")  




## Part C ##

## Table C1: models with ordinary least squares estimates
Model1b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN , data= data_long2, weights = wFactor)
Model2b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp , data= data_long2, weights = wFactor)
Model3b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + newer_party , data= data_long2, weights = wFactor)
Model4b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_LR  , data= data_long2, weights = wFactor)
Model5b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + DemBeh_01 , data= data_long2, weights = wFactor)
Model6b<-lm(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + newer_party + IdeolModeration_LR + DemBeh_01, data= data_long2, weights = wFactor)

tab_model (Model1b, Model2b, Model3b, Model4b, Model5b, Model6b, show.ci = FALSE, show.se = TRUE, digits.p = 2, emph.p =FALSE, file ="Table_C1.doc")


## Table C2: socio-political correlates of party legitimacy perceptions
Model1<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
tab_model (Model1, show.ci = FALSE, show.se = TRUE, emph.p =FALSE, file = "Table_C2.doc")


## Table C3: testing H3 with GAL-TAN scale
Model4a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_GT  +  (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model4b<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_LR + IdeolModeration_GT  +  (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)

tab_model (Model4a, Model4b, show.ci = FALSE, show.se = TRUE, emph.p =FALSE, digits.p = 2, file ="Table_C3.doc")


## Table C4: institutional trust as a moderator

# evaluating the index of trust
Trust<-data_long2[c("Q13_101", "Q13_102", "Q13_103", "Q13_104")]
cor(Trust, use = "pairwise.complete.obs") 

cronbach.alpha (Trust, na.rm=TRUE)

Trust_NA<-na.omit(Trust)
Trust_NA.pca<-prcomp(Trust_NA, center = TRUE, scale =TRUE)
summary(Trust_NA.pca)
Trust_NA.pca$rotation 

# produce input for Table C4
Model_T1<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Trust_01 + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model_T2<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + Trust_01 + Trust_01 * gov_exp + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model_T3<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + newer_party + Trust_01 + Trust_01 * newer_party + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model_T4<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_LR  + Trust_01 + Trust_01 * IdeolModeration_LR + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)
Model_T5<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + DemBeh_01 + Trust_01 + Trust_01 * DemBeh_01 + (1|Party) + (1|COUNTRY), data=data_long2, weights = wFactor)

tab_model (Model_T1, Model_T2, Model_T3, Model_T4, Model_T5, show.ci = FALSE, show.se = TRUE, digits.p = 2, emph.p =FALSE, file ="Table_C4.doc")

p1<-plot_model(Model_T2, type = "int", title = "", axis.title = c("governing experience", "legitimacy perceptions (0-1)")) + ggplot2::labs(colour = "Trust")
p2<-plot_model(Model_T3, type = "int", title = "", axis.title = c("party age", "legitimacy perceptions (0-1)")) + ggplot2::labs(colour = "Trust")
p3<-plot_model(Model_T4, type = "int",  title = "", axis.title = c("ideological moderation (0-1)", "legitimacy perceptions (0-1)")) + ggplot2::labs(colour = "Trust")
p4<-plot_model(Model_T5, type = "int",  title = "", axis.title = c("democratic behavior (0-1)", "legitimacy perceptions (0-1)")) + ggplot2::labs(colour = "Trust")

# all in one plot
jpeg("FigureC1.jpeg", width= 10, height= 8, units='in', res=300)
ggarrange(p1, p2, p3, p4, labels = c("A","B","C","D"), ncol = 2, nrow = 2)
dev.off()



## Table C5: controlling for party sympathy
Model1c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)
Model2c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + gov_exp + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)
Model3c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + newer_party + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)
Model4c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + IdeolModeration_LR  +  (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)
Model5c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + DemBeh_01 + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)
Model6c<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + Sympathy + gov_exp + newer_party + IdeolModeration_LR + DemBeh_01 + (1|Party) + (1|COUNTRY), data= data_long2, weights = wFactor)

tab_model (Model1c, Model2c, Model3c, Model4c, Model5c, Model6c, show.ci = FALSE, show.se = TRUE,digits.p = 2, emph.p =FALSE, file ="Table_C5.doc")



## Table C6: same sample throughout
data_fullM<-data_long2[c("Leg_01", "female", "Age", "Edu_o", "Extremity_SelfLRN", "gov_exp", "newer_party", "IdeolModeration_LR", "DemBeh_01", "Party", "COUNTRY", "wFactor")]
data_fullM<-na.omit(data_fullM)

Model1a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)
Model2a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)
Model3a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + newer_party + (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)
Model4a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + IdeolModeration_LR  +  (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)
Model5a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + DemBeh_01 + (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)
Model6a<-lmer(Leg_01~ female + Age + Edu_o + Extremity_SelfLRN + gov_exp + newer_party + IdeolModeration_LR + DemBeh_01 + (1|Party) + (1|COUNTRY), data= data_fullM, weights = wFactor)

tab_model (Model1a, Model2a, Model3a, Model4a, Model5a, Model6a, show.ci = FALSE, show.se = TRUE, digits.p = 2, emph.p =FALSE, file ="Table_C6.doc")


